C THIS PROGRAM IS DESIGNED FOR THE CALCULATION OF A GEOID UNDULATION
C AT A POINT WHOSE LATITUDE AND LONGITUDE IS SPECIFIED. THE PROGRAM
C IS DESIGNED TO USE THE POTENTIAL COEFFICIENT MODEL EGM96 AND A
C SET OF SPHERICAL HARMONIC COEFFICIENTS OF A CORRECTION TERM.
C THE CORRECTION TERM IS COMPOSED OF SEVERAL DIFFERENT COMPONENTS
C THE PRIMARY ONE BEING THE CONVERSION OF A HEIGHT ANOMALY TO A GEOID
C UNDULATION. THE PRINCIPLES OF THIS PROCEDURE WERE INITIALLY
C DESCRIBED IN THE PAPER: USE OF POTENTIAL COEFFICIENT MODELS FOR GEOID
C UNDULATION DETERMINATION USING A SPHERICAL HARMONIC REPRESENTATION
C OF THE HEIGHT ANOMALY/GEOID UNDULATION DIFFERENCE BY R.H. RAPP,
C JOURNAL OF GEODESY, 1996.
C THIS PROGRAM IS DESIGNED TO BE USED WITH THE CONSTANTS OF EGM96
C AND THOSE OF THE WGS84(G873) SYSTEM. THE UNDULATION WILL REFER TO
C THE WGS84 ELLIPSOID.
C SPECIFIC DETAILS ON THE UNDULATION COMPUTATION WILL BE FOUND IN THE
C JOINT PROJECT REPORT DESCRIBING THE DEVELOPMENT OF EGM96.
C THIS PROGRAM IS A MODIFICATION OF THE PROGRAM DESCRIBED IN THE
C FOLLOWING REPORT:
C A FORTRAN PROGRAM FOR THE COMPUTATION OF GRAVIMETRIC QUANTITIES FROM
C HIGH DEGREE SPHERICAL HARMONIC EXPANSIONS, RICHARD H. RAPP,
C REPORT 334, DEPARTMENT OF GEODETIC SCIENCE AND SURVEYING, THE OHIO
C STATE UNIVERSITY, COLUMBUS, 1982
C THIS PROGRAM WAS PUT IN THIS FORM IN DEC 1996.
C RHRAPP.F477.NONLY
C 
C
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  The input files consist of:
C
C                correction coefficient set ("CORRCOEF") => UNIT = 1
C                    potential coefficient set ("EGM96") => UNIT = 12
C                points at which to compute (INPUT.DAT") => UNIT = 14
C
C  The output file is:
C
C                    computed geoid heights ("OUTF477") => UNIT = 20
C
C    FILE ASSIGNMENT REVISIONS AT NIMA, DECEMBER 1996.
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C   DIMENSIONS OF P,Q,HC,HS MUST BE AT LEAST ((MAXN+1)*(MAXN+2))/2,
C   DIMENSIONS OF SINML,COSML,SCRAP MUST BE AT LEAST MAXN,
C        WHERE MAXN IS MAXIMUM ORDER OF COMPUTATION
C THE CURRENT DIMENSIONS ARE SET FOR A MAXIMUM DEGREE OF 360
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8  P(65341),SCRAP(361),RLEG(361),DLEG(361),RLNN(361),
     *SINML(361),COSML(361)
      REAL*8 HC(65341),HS(65341),CC(65341),CS(65341)
      DATA RAD/57.29577951308232D0/,HT/0.0D0/
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCC                                                           CCCCCCCCCCC
CCCCCCCCCC                 INPUT/OUPUT FILE NAMES                    CCCCCCCCCCC
CCCCCCCCCC                                                           CCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C CORRECTION COEFFICIENT FILE
      open(1,file='CORRCOEF',status='old')
C
C INPUT DATA FILE
      open(14,file='INPUT.DAT',status='old')
C
C  POTENTIAL COEFFICIENT FILE
      open(12,file='EGM96',status='old')
C
C  OUPUT FILE
      open(20,file='OUTF477.DAT',status='unknown')
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
cccc      CALL ERRSET(208,256,-1,2,IX,0)
cccc010   READ(5,900) NMAX
      NMAX = 360
900   FORMAT(I3)
      WRITE(6,910) NMAX
910   FORMAT(1H1,//' MAXIMUM DEGREE =',I4//)
      L=65341
      DO 71 I=1,L
      CC(I)=0.0D0
   71 CS(I)=0.0D0
C THE CORRECTION COEFFICIENTS ARE NOW READ IN
   72 READ(1,*,END=79)N,M,T1,T2
      IG=(N*(N+1))/2+M+1
      CC(IG)=T1
      CS(IG)=T2
      GO TO 72
   79 CONTINUE
C THE POTENTIAL COEFFICIENTS ARE NOW READ IN AND THE REFERENCE
C EVEN DEGREE ZONAL HARMONIC COEFFICIENTS REMOVED TO DEGREE 6
      CALL DHCSIN(NMAX,F,RJ2,RJ4,RJ6,HC,HS)
C SETTING IFLAG=1 PREVENTS LEGENDRE FUNCTION DERIVATIVES BEING TAKEN
C IN SUBROUTINE LEGFDN
      IFLAG=1
      IR=0
      K=NMAX+1
      FLATL=91.0D0
C READ GEODETIC LATITUDE,LONGITUDE AT POINT UNDULATION IS WANTED
cccc  030 READ(5,100,END=090) FLAT,FLON
  030 READ(14,*,END=090) FLAT,FLON
  100 FORMAT(6X,2F10.5)
C COMPUTE THE GEOCENTRIC LATITUDE,GEOCENTRIC RADIUS,NORMAL GRAVITY
      CALL RADGRA(FLAT,FLON,HT,RLAT,GR,RE)
      IF(FLATL.EQ.FLAT) GO TO 040
      RLAT1=RLAT
      RLAT=1.5707963267948966D0-RLAT
      FLATL=FLAT
      DO 25 J=1,K
      M=J-1
      CALL LEGFDN(M,RLAT,RLEG,DLEG,NMAX,IR,RLNN,IFLAG)
      DO 26 I=J,K
      N=I-1
      LOC=(N*(N+1))/2+M+1
  26  P(LOC)=RLEG(I)
  25  CONTINUE
040   RLON=FLON/RAD
      CALL DSCML (RLON,NMAX,SINML,COSML)
      CALL HUNDU(U,NMAX,P,HC,HS,SINML,COSML,
     *GR,RE,RLAT1,CC,CS,HACO)
C  U IS THE GEOID UNDULATION FROM THE EGM96 POTENTIAL COEFFICIENT MODEL
C  INCLUDING THE HEIGHT ANOMALY TO GEOID UNDULATION CORRECTION TERM
C  AND A CORRECTION TERM TO HAVE THE UNDULATIONS REFER TO THE
C  WGS84 ELLIPSOID. THE GEOID UNDULATION UNIT IS METERS.
cccc     WRITE(6,101)FLAT,FLON,U
      WRITE(20,101) FLAT,FLON,U
 101  FORMAT(2F14.7,F10.3)
      GO TO 30
  830 CONTINUE
 90   STOP
      END
      SUBROUTINE HUNDU(UNDU,NMAX,P,HC,HS,
     *SINML,COSML,GR,RE,ANG,CC,CS,HACO)
      IMPLICIT REAL*8 (A-H,O-Z)
cccc      DIMENSION SINML(1),COSML(1),P(1)
cccc      REAL*8 HC(1),HS(1),CC(1),CS(1)
      DIMENSION SINML(361),COSML(361),P(65341)
      REAL*8 HC(65341),HS(65341),CC(65341),CS(65341)
C  CONSTANTS FOR WGS84(G873)
      DATA GM/.3986004418D15/,AE/6378137.0D0/
C  GM IN UNITS OF M**3/S**2
      AR=AE/RE
      ARN=AR
      AC=0.0
      A=0.0
      B=0.0
      K=3
      DO 030 N=2,NMAX
      ARN=ARN*AR
      K=K+1
      SUM=P(K)*HC(K)
      SUMC=P(K)*CC(K)
      SUM2=0.0
      DO 020 M=1,N
      K=K+1
      TEMPC=CC(K)*COSML(M)+CS(K)*SINML(M)
      TEMP=HC(K)*COSML(M)+HS(K)*SINML(M)
      SUMC=SUMC+P(K)*TEMPC
020   SUM=SUM+     P(K) *TEMP
      AC=AC+SUMC
 30   A=A+SUM*ARN
      AC=AC+CC(1)+P(2)*CC(2)+P(3)*(CC(3)*COSML(1)+CS(3)*SINML(1))
      HACO=AC/100.D0
      UNDU=A*GM/(GR*RE)
C ADD HACO TO CONVERT HEIGHT ANOMALY ON THE ELLIPSOID TO THE UNDULATION
C ADD -0.53M TO MAKE UNDULATION REFER TO THE WGS84 ELLIPSOID.
      UNDU=UNDU+HACO-0.53D0
      RETURN
      END
      SUBROUTINE  DSCML  (RLON,NMAX,SINML,COSML)
      IMPLICIT REAL*8 (A-H,O-Z)
cccc      DIMENSION SINML(1),COSML(1)
      DIMENSION SINML(361),COSML(361)
      A=DSIN(RLON)
      B=DCOS(RLON)
      SINML(1)=A
      COSML(1)=B
      SINML(2)=2.0*B*A
      COSML(2)=2.0*B*B-1.0
      DO 010 M=3,NMAX
      SINML(M)=2.0*B*SINML(M-1)-SINML(M-2)
010   COSML(M)=2.0*B*COSML(M-1)-COSML(M-2)
      RETURN
      END
      SUBROUTINE DHCSIN (NMAX,F,J2,J4,J6,HC,HS)
      IMPLICIT REAL*8 (A-H,O-Z)
      CHARACTER IC
      REAL*8 J2,J4,J6,J8,J10,K
      REAL*8 HC(65341),HS(65341)
C  CONSTANTS FROM WGS84(G873)
      DATA GM,OMEGA,AE/.3986004418D20,7.292115D-5,6378137.0D0/
      DATA RF/298.257223563D0/
C THE EVEN DEGREE ZONAL COEFFICIENTS GIVEN BELOW WERE COMPUTED FOR THE
C WGS84(G873) SYSTEM OF CONSTANTS AND ARE IDENTICAL TO THOSE VALUES
C USED IN THE NIMA GRIDDING PROCEDURE. COMPUTED USING SUBROUTINE
C GRS WRITTEN BY N.K. PAVLIS
      J2=0.108262982131D-2
      J4=-.237091120053D-05
      J6=0.608346498882D-8
      J8=-0.142681087920D-10
      J10=0.121439275882D-13
      M=((NMAX+1)*(NMAX+2))/2
      DO 001 N=1,M
      HC(N)=0.0
001   HS(N)=0.0
  912 FORMAT(2I4,2D30.20)
cccc 02   READ(12,912,END=3)N,M,C,S
 02   READ(12,*,END=3)N,M,C,S,EC,ES
      IF(N.GT.NMAX) GO TO 002
      N=(N*(N+1))/2+M+1
      HC(N)=C
      HS(N)=S
      GO TO 002
3     HC(4)=HC(4)+J2/DSQRT(5.D0)
      HC(11)=HC(11)+J4/3.0D0
      HC(22)=HC(22)+J6/DSQRT(13.D0)
      HC(37)=HC(37)+J8/DSQRT(17.D0)
      HC(56)=HC(56)+J10/DSQRT(21.D0)
      RETURN
      END
      SUBROUTINE LEGFDN(M,THETA,RLEG,DLEG,NMX,IR,RLNN,IFLAG)
C
C            THIS SUBROUTINE COMPUTES  ALL NORMALIZED LEGENDRE FUNCTION
C            IN "RLEG" AND THEIR DERIVATIVES IN "DLEG". ORDER IS ALWAYS
C            M , AND COLATITUDE IS ALWAYS THETA  (RADIANS). MAXIMUM DEG
C            IS  NMX  . ALL CALCULATIONS IN DOUBLE PRECISION.
C            IR  MUST BE SET TO ZERO BEFORE THE FIRST CALL TO THIS SUB.
C            THE DIMENSIONS OF ARRAYS  RLEG, DLEG, AND RLNN  MUST BE
C            AT LEAST EQUAL TO  NMX+1  .
C
C            THIS PROGRAM DOES NOT COMPUTE DERIVATIVES AT THE POLES .
C
C            IF    IFLAG = 1  , ONLY THE LEGENDRE FUNCTIONS ARE
C            COMPUTED.
C
C      ORIGINAL PROGRAMMER :OSCAR L. COLOMBO, DEPT. OF GEODETIC SCIENCE
C      THE OHIO STATE UNIVERSITY, AUGUST 1980 . ******************
C
C
      IMPLICIT REAL*8 (A-H,O-Z)
cccc      DIMENSION RLEG(1),DLEG(1),RLNN(1)
      DIMENSION RLEG(361),DLEG(361),RLNN(361)
     2, DRTS(1300),DIRT(1300)
      NMX1 = NMX+1
      NMX2P = 2*NMX+1
      M1 = M+1
      M2 = M+2
      M3 = M+3
      IF(IR.EQ.1) GO TO 10
      IR = 1
      DO 5     N = 1,NMX2P
      DRTS(N) = DSQRT(N*1.D0)
    5 DIRT(N) = 1.D0/DRTS(N)
   10 COTHET = DCOS(THETA)
      SITHET = DSIN(THETA)
      IF(IFLAG.NE.1.AND.THETA.NE.0.D0)SITHI = 1.D0/SITHET
C
C            COMPUTE THE LEGENDRE FUNCTIONS .
C
      RLNN(1) = 1.D0
      RLNN(2) = SITHET*DRTS(3)
      DO 15    N1 = 3,M1
      N = N1-1
      N2 = 2*N
   15 RLNN(N1) = DRTS(N2+1)*DIRT(N2)*SITHET*RLNN(N1-1)
      IF(M.GT.1) GO TO 20
      IF(M.EQ.0) GO TO 16
      RLEG(2) = RLNN(2)
      RLEG(3) = DRTS(5)*COTHET*RLEG(2)
      GO TO 20
   16 RLEG(1) = 1.D0
      RLEG(2) = COTHET*DRTS(3)
   20 CONTINUE
      RLEG(M1) = RLNN(M1)
        IF(M2.GT.NMX1) GO TO 35
      RLEG(M2) = DRTS(M1*2+1)*COTHET*RLEG(M1)
        IF(M3.GT.NMX1) GO TO 35
      DO 30     N1 = M3,NMX1
      N = N1-1
      IF(M.EQ.0.AND.N.LT.2.OR.M.EQ.1.AND.N.LT.3) GO TO 30
      N2 = 2*N
      RLEG(N1) = DRTS(N2+1)*DIRT(N+M)*DIRT(N-M)*(DRTS(N2-1)*COTHET*
     2 RLEG(N1-1)-DRTS(N+M-1)*DRTS(N-M-1)*DIRT(N2-3)*RLEG(N1-2))
      GO TO 30
   30 CONTINUE
   35 IF(IFLAG.EQ.1) RETURN
      IF(SITHET.EQ.0.D0) WRITE(6,99)
   99 FORMAT(//' *** LEGFDN  DOES NOT COMPUTE DERIVATIVES AT THE POLES
     2 *****************'//)
      IF(SITHET.EQ.0.D0) RETURN
C
C            COMPUTE ALL THE DERIVATIVES OF THE LEGENDRE FUNCTIONS.
C
      RLNN(1) = 0.D0
      RLN = RLNN(2)
      RLNN(2) = DRTS(3)*COTHET
      DO 40      N1 = 3, M1
      N = N1-1
      N2 = 2*N
      RLN1 = RLNN(N1)
      RLNN(N1) = DRTS(N2+1)*DIRT(N2)*(SITHET*RLNN(N)+COTHET*RLN)
      RLN = RLN1
   40 CONTINUE
      DLEG(M1) = RLNN(M1)
        IF(M2.GT.NMX1) RETURN
      DO 60      N1 = M2,NMX1
      N = N1-1
      N2 = N*2
      DLEG(N1) = SITHI*( N   *RLEG(N1)*COTHET-DRTS(N-M)*DRTS(N+M)*
     2 DRTS(N2+1)*DIRT(N2-1)*RLEG(N))
   60 CONTINUE
      RETURN
      END
      SUBROUTINE RADGRA(FLAT,FLON,HT,RLAT,GR,RE)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 N,K
C THIS SUBROUTINE COMPUTES GEOCENTRIC DISTANCE TO THE POINT,
C THE GEOCENTRIC LATITUDE,AND
C AN APPROXIMATE VALUE OF NORMAL GRAVITY AT THE POINT BASED
C THE CONSTANTS OF THE WGS84(G873) SYSTEM ARE USED
      DATA AE/6378137.D0/,E2/.00669437999013D0/,RAD/57.29577951308232D0/
      DATA GEQT,K/9.7803253359D0,.00193185265246D0/
      FLATR=FLAT/RAD
      FLONR=FLON/RAD
      T1=DSIN(FLATR)**2
      N=AE/DSQRT(1.-E2*T1)
      T2=(N+HT)*DCOS(FLATR)
      X=T2*DCOS(FLONR)
      Y=T2*DSIN(FLONR)
      Z=(N*(1.-E2)+HT)*DSIN(FLATR)
      N=AE/DSQRT(1.-E2*T1)
C COMPUTE THE GEOCENTRIC RADIUS
      RE=DSQRT(X**2+Y**2+Z**2)
C COMPUTE THE GEOCENTRIC LATITUDE
      RLAT=DATAN(Z/DSQRT(X**2+Y**2))
C COMPUTE NORMAL GRAVITY:UNITS ARE M/SEC**2
      GR=GEQT*(1.+K*T1)/DSQRT(1.-E2*T1)
      RETURN
      END
